#############################################################################################
#############################################################################################
####Non-state Atrocities in Capital Cities - Insurgent Atrocities (1996-2009) Global Map ####
#############################################################################################
#############################################################################################
library(foreign)
library(doBy)
library(plyr)
library(ggplot2)
library(ggthemes)
library(maps)
library(maptools)
library(ggmap)
library(mapdata)
library(RgoogleMaps)

# Set working library
setwd("~/Data/Global Analysis/")
# Read in main data
main.data <- read.dta("full.grid.upd.dta")

###Create data for plotting
#Sum insurgent atrocities at year t by coordinates
values<- summaryBy(incidentnonstatefull ~ latitude + longitude, keep.names = T, data=main.data, FUN=c("sum"), na.rm=TRUE)
#Create logged values
values$incidentnonstatefull <- log(values$incidentnonstatefull+1)
#Normalize values by relations to the grid cell with the highest number of atrocities in the sample
values$incidentnonstatefull <-values$incidentnonstatefull/max(values$incidentnonstatefull)
#Create a color spectrum for the map
col.s <- rgb(1, 0, 0, alpha=1*values$incidentnonstatefull, maxColorValue = 1)
#Create color spectrum for the legend
c0 <- rgb(1, 1, 1, maxColorValue = 1)
c1 <- rgb(1, 0, 0, alpha=0.25, maxColorValue = 1)
c2 <- rgb(1, 0, 0, alpha=0.50, maxColorValue = 1)
c3 <- rgb(1, 0, 0, alpha=0.70, maxColorValue = 1)
colvec.at <- rbind(c0, c1, c2, c3)
#Create categories for the legend
leg.at.txt <- c("No Atrocities", "Low Frequency", "Medium Frequency", "High Frequency")
###Plot the map
#pdf("atrocitymap.pdf")
map("worldHires", col="white", fill=TRUE,xlim=c(-180, 180), ylim=c(-58, 90))
points(values$longitude, values$latitude, pch=16, col=col.s, cex=0.7)
legend("bottomleft",title = "Variation in The Frequency of Atrocities", 
       legend = leg.at.txt, 
       fill=colvec.at,
       cex = 0.56,
       bty = "n")
#dev.off()